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AN ALGORITHM FOR COMPUTING SPHERICAL PARTIALS 

Raymond V Borchers 
Program Systems Branch 
Mission & Trajectory Analysis Division 

ABSTRACT 

This paper presents a rapid method for computing satellite 
accelerations, position partials, and partials with respect to har- 
monic coefficients in the earth’s geopotential using spherical re- 
currences. Some timing estimates and accuracy comparisons are 
given as a function of order of the harmonics for two different 
satellites. 
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AN ALGORITHM FOR COMPUTING SPHERICAL PARTIALS 

\ , 

INTRODUCTION 

The nonlinear equations of motion of a satellite moving in the gravitational 
field of the earth are presently being solved by high order predictor-corrector 
methods of numerical integration. The advent of extremely accurate electronic 
and optical measurements makes it both desirable and necessary to include in 
the force model terms of high degree and order as well as other forces such as 
drag, solar radiation, etc. when estimating geodetic parameters, tracking station 
coordinates, etc. 

This paper describes a fast method for computing accelerations, partial 
derivatives with respect to position, and partial derivatives with respect to 
harmonic coefficients in the earth’s geopotential using recurrences in spherical 
coordinates. This method (1], [3] was compared with De Witt’s method [ 2 J 
which uses recurrences in cartesian coordinates. The spherical recursion 
algorithm was found to be as accurate and at least twice as fast in t 1 <3 calculation 

t 

of accelerations alone (1]. The author reached the same conclusion and, in ad- 
dition found the spherica 1 version much simpler from both analysis and program- 
ming points of view. The method described here computes the total acceleration 
(Keplerian plus perturbative) due to the geopotential field for an artificial satel- 
lite orbiting the earth. The principal features of this method have been applied 

i 

in a program which estimates geodetic parameters, tracking station coordinates, 
etc. using multiple arcs of various satellites. 
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GRAVITATIONAL POTENTIAL OF THE EARTH 

The potential function of the earth (adopted at the meeting of Commission 7 
on celestial mechanics at the Berkeley meeting of the International Astronomical 
Union in August 1961) is given by 



where, 

P n m (^) = associated Legendre functions of degree n and order m, 

M = z/r = sine of the geocentric latitude, 

G = Universal gravitational constant, 

M = mass of the earth, 

R = mean equatorial radius of the earth, 

\ = east longitude of the satellite, 
r = geocentric distance to satellite, 

C n m , S n m = unnormalized harmonic coefficients. 

If the center of gravity of the earth is chosen as the center of the coordinate 
system, the terms for n = 1 do not exist, i.e., c® = c t l = s ® = s t l = 0. The 
form of the potential then becomes 
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whore 


N maximum order of harmonic expansion selected. 

It is well known that the potential function U is a solution of Laplace's 
equation, which in rectangular coordinates is given by 


a 2 u a 2 u a 2 u 

ax 2 ay 2 a z 2 


and in spherical coordinates by 


a 2 u 1 <? 2 u i <?u i «j 2 u 

av 1 r *r r 2 d0 2 r 2 tanV '<^ r 2 cos 2 «// «9\ 2 = °* 


( 2 ) 


(3) 


COOKDINATE SYSTEMS 

Two different orthogonal right handed coordinate systems shown in Figure 1 
are used. These systems are 

(1) An inertial coordinate system referenced with respect to the first point 
of Aries, T, and 

<2) A geocentric coordinate system referenced with respect to the Green- 
wich meridian, 

where 


r 

T 


x. y. z 
X G ’ ^ G ’ Z G 


satellite position ve :cor in either coordinate system, 
first point of Aries, 

inertial cartesian coordinates of position vector r, 
geocentric coordinates referenced to Greenwich meridian. 
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0 


z, z c 



Figure 1. 


a - right ascension of the satellite, 

» GM = right ascension of the Greenwich meridian, 
k -■ east longitude of the satellite, 

0 = giocentric latitude of the satellite. 

We obtain the angle k and the sine of 0 from the following equations 

y z 

a " toi~ l — , k = a - a Q , s in 0 ~ “ * 

The partials of r, 0, k with respect to x, y, z can be obtained from differ- 
entiating the following expressions: 


2 - 


; 2 + y 2 + Z 2 


0 = tan 


- 1 


(x 2 + y 2 ) 


1/2 


K = - 


a Q + tan 


-l JL 


f 
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We obtain after differentiation the following partial derivatives: 


d r X 
r • 


dip ~ xz 

'**’ r J (x J + y J ) l/3 ' 


«9\ ~ y 

<?x y 2 + y 2 


<)y r • 

dip - yz 

17 r JJ x l + y J)l/2 


X 



d T z_ 

7z r • 


6 Kp (x 2 ♦ y 2 '| 

77 = J — 72 — 


SK 

Tz = °* 


fhe equations for the total accelerations hen become 


x dU y 

£U 

xz f>U 

t dr " x 2 4 y 2 

d A r 2 

J x 2 4 y 2 jl 2 ft . 

y dU x 

* 9 U 

yz 

rTr x 2 4 y 2 

TK “ r 2 

( x 2 4 y 2 J 1/2 

z /x 2 + y 

2 <?u 


r ~d r r 2 

” 17 ' 



It should be noticed that in the equations for the partials of U with respect 
to the spherical coordinates r , 0, and K certain quantities such as P n m , C n m cos m\ 

+ S n m sin mX, and GM r are common to these partials and need be computed only 
once. Furthermore, recursive relations are used for generating expressions 
such as sin nv\, cos m\, m tan \p 9 and P n m so as to increase the speed of computation. 
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The formulas used In recursion are: 




P n °(h) = [( 2n - 1 ) sin <// P n °- , (m ) - (n - 1 ) P n °. 2 (m)] In 

P n m 00 = P n "- 2 + ( 2n - 1 ) cos \p P"-V (m) w 7 0, m < n 
For sectorials m n and we have 

P," (M) = (2n - 1) cos <// P n n .‘, ! (M)m ^ 0 , m = n . 

sin mV = '/ cos \ sin (m ~ 1 ) \ - sin (m * 2)\ 

cos m\ : 2 cos \ cos (m - 1 ) \ - cos (m - 2) \ 

mtani/r = £(m - l ) tan t/'J + tnrwp . 

ACCELERATIONS 

The total gravitational acceleration vector F 0 in inertial coordinates (x, y t z) 
is given by 

dU <?U d\J 
F 0 = l - r?x * 1 dy + ^ dz • 

where i, j , k form a triad of mutually orthogonal unit vectors in the inertial co- 
ordinate system, and the partial derivatives of the potential function U are given 
with respect to inertial coordinates x , y , z . Applying the chain rule of calculus 
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to the function 


U " U(x, y, z) = U(r, 0, \) , 

we obtain the following equations for the accelerations: 


x 


<HJ_ 

<hj 

a r 

on 

aj. 

<9U 

a x 

<9x 

a r 

ax 

* <90 

ax 

4 ax 

<9x ’ 

<9U 

<HJ 

a r 

gu 

a^p 

<9U 

dX 

dy 

a r 

a y 

4 <90 

ay 

4 <9X 

ay ’ 

an 

a v 

a r 



<9U 

ax 

a z 

a r 

a z 

* <90 

a z 

4 ax 

<9z * 


(4) 

The partials jf U with respect to r, «/ , > , i.e., <) U /<? r , <9U <90, au^ax are ob- 
tained by explicit differentiation of the potential function U and are expressed by 
the following equations: 


rHJ 
<9 r 


= - 7 (t 1 ) 1 + (t)” <" + «> £ ( C COS "* 4 S n” si " mA ) P ; 

n - 2 m= 0 


<90 = ( G r M ) 2^ (f) 22 ( C, ' m COS + S " Sin [ Pn, " + 1 ' m tan ^ P n"] 

n - 2 m * 0 

W = (t) F. (r) 2Z m ( S "” cos " C '’ m sin mX ) p ” ' 

n = 2 m ~ 0 


(5) 





POSITION PAPTIALS 
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The desired partial derivatives with respect to inertial cartesian coordinates 
can be written in the following matrix form 


a 2 u 

a 2 u 

a 2 u 


ax 

ax 

ax 

ax 2 

dx dy 

dx a z 



<?y 

dz 

a 2 u 

a 2 u 

a 2 u 


ay 

ay 

dy 

dy dx 

ay 2 

dy dz 



37 

71 

,? 2 U 

a 2 u 

a 2 u 


a*z 

az 

dz 


dzdy 

az 2 


J77 

57 

Jz_ 


From Laplace's Equation 2 we can solve for d 2 U/dx 2 , i.e., 

<9 2 U / a 2 U + d 2 U \ 

dx 2 \ a y 2 dz 2 ) 


Since the above matrix is symmetric, we need only compute the elements 
above the principal diagonal and two remaining elements on the principal diag- 
onal yielding a total of 5 elements. 


The elements in the above matrix can be obtained by differentiating the ex- 
pressions (4) with respect to x, y, and z and are given as follows: 


dx 

7x 


au aj_r 
dr ax 2 


au o 2 ip au d 2 k a 2 u /ar\ 2 a 2 u/a^\ 2 a 2 u /a\\ 2 

+ 


/ a 2 u ar dy a 2 u ar a\ 

+ ^ \d r dyp a x a x + d r d\ dx dx 


a 2 u a^ 3k\ 

+ dypdh dx dxj 
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t 


* 


7 x 7U 7 I 2 r 7U 7 2 v 7U 9 2 k 

7y 77 7y 7x ' 70 7y77 * 77 7Ty~dx 

7 2 U dr dr 7j_U 70 70 d 2 U dK dk 
dr 2 Jy 7x 4 ,902 7y 7x 4 ^2 7y 77 


7 2 U (dr d*p dr 7 0 \ 7 2 U Id r dK dr dk \ 

+ 77“70 l7P 77 + 77 T^j + 7777 \7£ 77 * 77 3y / 


7 2 u /70 7\ M>dk\ 

+ 7R7 V7y 77 4 57 7^ / * 


ax au a 2 r au a*^ au a*x 

77 77 dx dz * 77 7x7 z ' 7X 7777 


7 2 U7r7r 7 2 U 70 70 7 2 U 7\ 7\ 

<9 r 2 ^7 ' Tx * 4 70 2 fJz + 7\ 2 <* z 


7 2 U 

/ 7 r 70 

7 r 

70 \ 

7 2 U 

/ 7 )• 7A. 

7 r 

7A. \ 

7 r 70 

\77 77 

7x 

77) 

1 7777 

\7Y 77 

+ 77 

77/ 


7 2 U /70 7 A. 
+ 70 7 A. \ 7 z 7 x 


70 7\ \ 

+ 77 77/ 


7y 7U 7 2 r 

7y 77 dy 2 


7U 7 2 0 7U 7 2 \ 7 2 U f d r \ 2 7 2 U ^70 \ 2 7HJ ^7\\ 2 

5 y 2 ’ 3X 2 ‘ ^ r 2 UyJ + d4 ji \Xy> * aK 2 \7yl 

I 7 2 U 7 r 70 7 2 r 7r 7A. 7 2 U 70 7 A. \ 

^ v 77T0 7y 7y 4 7r 7 v 7y 7y 4 70 7 A. 7y 7y / 
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s 


<9y (9U <9 2 r d[) <9 2 p <9U <9 2 k 

dz Hr dy Hz * dy Hz * 37" dy dz 


d 2 \] dr dr <? 2 U dy dp d 2 U dK dk 
dr 2 ' dip 2 ' d\ 2 ** 


(9 2 U (dr dip dr dp\ <9 2 U (dr dk dr dk\ 
+ dTTJ, 37 ' Ti 'Sy) 4 37^A \TF,Jl'Tz^) 


d 2 U (d4> dk dip dk\ 
4 £p dk \dy d z + d z dy) 


dz dV d 2 r d\] d 2 ip dV d 2 k d 2 U (dr \ 2 d 2 [) (dip\ 2 d 2 U (dk\ 2 

^ ^ dz 2 *5 dz 2 ^ dz 2 ~d7 2 + Ip 2 + *\ 2 


( d 2 U <9r dip d 2 U ar <9\ (9 2 U (90 <9\\ 

4 “ V/9 r (9v dz dz + d r d\ dz 37 4 (9 0 (9 A 3z 37 / 


We need the expression for the double gradient of U given by 



[jL] 

dr 




(9 2 U 
<9r 2 

(9 2 U 
dr dp 

(9 2 U 
d r dk 

7(VU) 

_d_ 

dp 

r (9u 
Ur 

<9U 

dp 

ii 

r i 

fC |<T> 

d 2 U 
dp d r 

d 2 U 
dp 2 

<9 2 U 
dp dk 


d 

Jk_ 




<9 2 U 
dk d r 

(9 2 U 
<9\ dp 

d 2 U 
dk 2 __ 


We notice that the 3x3 matrix is symmetric and hence need only to compute 
the elements on and above the principal diagonal. But we can also obtain one of 
the elements on the principal diagonal from Laplace’s Equation 3. This means 
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that only a total of 5 partials are requirod In this matrix. The second partials 
of U with respect to r, ip 9 and k are obtained by differentiating the expressions (5) 
for the first partials and are given by 


<9 2 U 
<9 r 2 


<9 2 U 
<9 r dip 


<9 2 U 
d r <5 a . 


(9 2 U 
34j 2 


d 2 {] 

dpi dk 


<9 2 u 

Ik 2 


CM 

.3 


N n n 

2 + ^ ( 7 ) (n + 2) (n ♦ 1) ( C n m cos mk + S n w sinm^P^ 

_ n :, 2 m= n 


N n n 

= 1(n+i) 2 ^( c "” cos '^ ts ” sin '^)( p »"* , ' mtan ' / ' p "") 
n=2 m=0 


N n 

GM r / r \ n i 1 

- Js 2_.' 7 ' (ntl) / ^(S’cosn* -C^simi*)?,,* 

n= 2 1 ~n 


N n 


= (t^)£](t) 7^ (c" cosh* 


n r 2 m~ 0 


+ S n m sin nv\) |tan 0P n m + 1 + [ m sec2 0 " m tan2 0 ” n ( n + 1 )] P,," 1 } 


N n 


= f. m ( s n l " cos ~ C n m »inn*)(P n "* 1 ~ m tan ^P n m ) 


n= 2 m= 0 


= (~ T~) 7^ ( 7 ) 7"! m2 ( C n" COS * S " sinm ^) p : 


n- 2 m= 0 
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We mentioned above that we needed only 5 partials. Hence, we arbitrarily 


chose to computed 2 U/d</> 2 from Laplace's Equation 3, i.e., 


d 2 U 


dU 

tan ^ 


I f 2 <J 2 u 

5 T ^ 


+ 2r 77 + 


1 d 2 il l 

i 2 0 d\ 2 J 


The second partial derivatives of r, and A with respect to x, y, and z can 
be represented in the symmetric matrix form 


d 2 r 


d 2 r 
dxdy 


d 2 r 
dx d z 


d 2 r 
dy dx 


d 2 r 


d 2 r 
dy d z 


d 2 r 
dz dx 


d 2 r 
d z dy 


d 2 r 


where the elements of the r matrix are given by 


d 2 r 


d 2 r 


d 2 r 


d 2 r 


d 2 r _ d 2 r 
dx dz d z dx 


d 2 r 


d 2 r 


1 


/ 


The partial derivatives of 0 with respect to x, y, and z can be represented by 
the symmetric matrix 


d 2 4> 

<9 2 0 

<9 2 0 

dx 2 

dxdy 

<9x <?z 

d 2 0 

d 2 0 

<9 2 0 

f9y <9x 


<9y <9z 

<9 2 0 

<9 2 0 

<9 2 0 

dz <9x 

dzdy 

5 z 2 


where the elements of the 0 matrix are given by 


<9 2 0 
dx 2 


4 (x 2 + y 2 ) 3 2 


[2x J ( 


) - r J y : 


<9 2 0 < 9 2 0 

dx dy d y (lx 


xyz 

r 4 (x 2 + y 2 ) 3 ' 2 


[j3(x 2 * y 2 ) * z 2 j 


<9 2 (9 2 0 

dx dz d zdx 


x 

r 4 (x 2 * y 2 ) l/2 


[z 2 -(x 2 +y 2 )] 


<9 2 0 
(9y 2 


z 

r 4 (x 2 + y 2 ) 3 2 




x 


*] 


<9 2 0 
c)y <9 z 


< 9 2 0 

<9 z 2 


jdy = r -4 (x 4 y2 )i72 [* 2 - (x 2 + y 2 )] 

- ^ ( x 2 + y 2 ) 1/2 . 
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The partial derivatives of k with respect to x, y , and z can be represented 
by the symmetric matrix 


d 2 k 

dx 2 

d 2 k 
(ty dx 

d 2 k 

d Z dx 


where the elements of the k matrix 


d 2 k d 2 k 

2xy 

dx 2 <9y 2 

(x 2 4y!)! • 

*9 2 \ 

d 2 k d 

d 7 2 

r 9 X d Z 03 


d 2 k 
dx dy 

d 2 k 

dxdz 


d 2 k 

^y 2 " 

d 2 \ 

dy dz 

t 

d 2 A, 
dz d y 

d 2 k 

dz 2 


are given by 


d 2 k 

d 2 k 

1 

dx dy 

dydx 

( x 2 4 y 2 ) 2 

ii 

& 

d 2 k 
dy dz 

d 2 k 

TTTiy 0 • 



The higher accuracy that is required in computing accelerations is not 
needed in the calculation of the position partials so that fewer terms are in- 
cluded in the geopotential. 


HARMONIC COEFFICIENT PARTIALS 

The equations for computing partial derivatives of acceleration with respect 
to harmonic coefficients C n m or S n m obtained by explicit differentiation of the 
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equations for the accelerations (4) are given by 


where 
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1 /GM\ / R\ n 

“7Y7"Mt) (n + l)cosmAP; 
(^t)(t) cosmX(P n mM -mtan^P") 


(9 

r9C m 

n 


<9 

r9S m 

n 



/GM\/R\ n 

“ \Tl\T) ^sinm\P n m 

1/GM\/R\ n . , 

7 \“r~ / \T/ ( n + 1) sinmXP; 
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s 



(?) (t) sinraAlp ;* 1 - m tan V-P n ") 


(?)(t) m«-o S m\P n " 


TIMING ESTIMATES 

Several runs were made on the IBM 360 model 91 which uses the OS MVT 
multiprogramming system. Some estimates of running time were obtained by a 
pair of routines called ETIMIN, ETIMOT (TIMEIN, TIMOUT). The first routine 
reads the system clock but returns no output and the second routine returns the 
time elapsed since the first routine was called. This actual elapsed time may 
include other users’ time as well. Consequently, the estimates do vary some- 

i 

what and are usually greater than they should be. 

The program was called upon repeatedly for 10,000 cases to simulate the 
calculations required for the same number of time points t. Some timing estimates 
are presented for accelerations plus position partials, and position partials only 
in Tables 1 and 3. The number of ’’significant” digits indicated in Tables 2 and 
4 represents the number of digits agreeing with the standard case (degree = 15, 
order = 15). The first and second columns in each table show the order and de- 
gree of the harmonic coefficients respectively for the computation of accelera- 
tions. In addition the third column provides the degree of the expansion used in 
computing the position partials. 
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Table 1 

Time Estimates 


Order 

Degree 

(accelerations) 

Degree 

(position 

partials) 

Running Time in Seconds 

Accelerations 

Accelerations Plus 
Position Partials 

15 

15 

2 

57.0 

59.9 

15 

15 

4 

57.0 

60.9 

15 

15 

6 

57.0 

62.0 

15 

15 

8 

57.0 

63.3 

15 

15 

10 

57.1 

64.1 

15 

15 

12 

57.1 

64.7 

15 

15 

15 

57.2 

65.5 

10 

10 

10 

32.3 

38.2 

5 

5 

5 

11.7 

14.5 

2 

2 

2 

3.4 

4.7 


TETR-C DATA 

x = -0.197 171 190 (06) Meters 
y = -0.646 084 338 (07) Meters 
z = +0.250 067 586 (07) Meters 
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Table 2 


Significant Digits 


Order 

Degree 

(accelerations) 

Degree 

(position 

partials) 

No. of Significant Digits 

Accelerations 

Accelerations Plus 
Position Partials 

15 

15 

2 

16 

3 

15 

15 

4 

16 

3 

15 

15 

6 

16 

4 

15 

15 

8 

16 

4 

15 

15 

10 

16 

5 

15 

15 

12 

16 

5 

15 

15 

15 

16 

16 

10 

10 

10 

6 

5 

5 

5 

5 

4 

3 

2 

2 

2 

4 

3 


TETR-C DATA 

x = -0.197 171 190 (06) Meters 

y = -0.646 084 338 (07) Meters 

z = +0.250 067 586 (07) Meters 
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Table 3 

Time Estimates 


Order 

Degree 

(accelerations) 

Degree 

(position 

partials) 

Running Time in Seconds 

Accelerations 

Accelerations Plus 
Position Partials 

Position 

Partials 

15 

15 

2 

45.5 

48.4 

46.6 

15 

15 

4 

45.5 

48.9 

48.0 

15 

15 

6 

45.5 

51.5 

48.4 

15 

15 

8 

45.5 

54.5 

49.5 

15 

15 

10 

45.5 

56.3 

49.9 

15 

15 

12 

45.5 

55.3 

50.4 

15 

15 

15 

45.5 

53.1 

52.7 

10 

10 

10 

26.2 

32.8 

31.1 

5 

5 

5 

— 

12.9 

12.3 

2 

2 

2 

— 

4.6 

4.4 


GEOS DATA 

x = +0.569 053 863 879 2412 (07) Meters 

y = +0.147 453 452 873 1973 (07) Meters 

z = +0.601 344 521 360 5027 (07) Meters 
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Table 4 


Significant Digits 


Order 

Degree 

(accelerations) 

Degree 

(position 

partials) 

No. of Significant Digits 

Accelerations 

Accelerations Plus 
Position Partials 

15 

15 

2 

16 

4 

15 

15 

4 

13 

5 

15 

15 

6 

16 

5 

15 

15 

8 

16 

5 

15 

15 

10 

16 

5 

15 

15 

12 

16 

6 

15 

15 

15 

16 

16 

10 

10 

10 

6 

5 

5 

5 

5 

6 

5 

2 

2 

2 

5 

4 


GEOS DATA 

x = +0.569 053 863 879 2412 (07) Meters 

y = +0.147 413 452 873 1973 (07) Meters 

z = +0.601 344 521 360 5027 (07) Meters 
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Hence, it is possible to compute position partials of the same or lower order 
than the accelerations. This option is available when position partials need not 
be computed as accurately as the accelerations. 
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